#### Map the wells in the dataset ----
use.png = TRUE

map.file.name = '/Well_Map_'%&%today()
map.file.name = map.file.name%&%ifelse(use.png,'.png')

dt.map = dt[!(longitude==0 | latitude==0)]
locs.map = dt.map[, .(longitude, latitude)]

# Create SpatialPoints object with coords and CRS
points_sp <- SpatialPoints(coords = locs.map,
                           proj4string = CRS("+proj=longlat +datum=WGS84"))

dt.map[fed.land==1, fed.land.factor := 'Federal']
dt.map[fed.land==0, fed.land.factor := 'Nonfederal']
dt.map[, fed.land.factor := as.factor(fed.land.factor)]
dt.map[, fed.type := as.factor(prod_type%&%fed.land.factor)]

dt.map[, levels(fed.type)] # Colors need to be in this order
col_lookup = c('GasFederal'='#04273C', 'GasNonfederal'='#88C4F4',
               'OilFederal'='#FF6663', 'OilNonfedereal'='#EAD367')
col_lookup = alpha(col_lookup, 0.4)
names(col_lookup) = rep(c('Gas','Oil'), each=2)%&%c('Federal','Nonfederal')
col_lookup[dt.map[26879:26880,fed.type]] # check it gets the right col

pch_lookup = copy(col_lookup)
class(pch_lookup) = 'numeric'
pch_lookup[1:2] = 4 # cross
pch_lookup[3:4] = 20 # circle

pch.cex_lookup = copy(pch_lookup)
pch.cex_lookup[1:2] = 0.7
pch.cex_lookup[3:4] = 0.9


# Create SpatialPointsDataFrame object
points_spdf <- SpatialPointsDataFrame(coords = locs.map,
                                      data = dt.map[,.(latitude, longitude, 
                                                       drill_type=as.factor(drill_type), 
                                                       prod_type=as.factor(prod_type),
                                                       fed.land.factor,
                                                       fed.type,
                                                       boe_prac_ip)],  
                                      proj4string = CRS("+proj=longlat +datum=WGS84"))

vessel = st_read(data.dir%&%'/Maps/fedlandp.gdb')
vessel = as(vessel,"Spatial") # convert to spatial point format

# You may need to run the following 1-2 lines if the ne_states() function call fails.
# install.packages('devtools')
# devtools::install_github("ropensci/rnaturalearthhires")
states_sp <- ne_states()

states_sp@proj4string # uses WGS84 datum, same as above. important

# dev.off()
if (!use.png) pdf(output.dir%&%map.file.name, family='serif', width=11, height=8.5)
if (use.png) {
  png(output.dir%&%map.file.name, family='serif', 
      width=12*1.3, height=12, units='in',
      res=500, type='cairo')
}
# Set margins for bottom, left, top, and right of plot
par(mar = c(1, 1, 3, 1))

# states_sp@bbox <- bbox(points_spdf) # include Alaska
states_sp@bbox <- matrix(c(-125, 26, -65, 49), nrow=2) # cut Alaska

plot(states_sp, 
     lwd=ifelse(use.png, 1.2, 1),
     col = '#FFFFFF',
     border = '#687D8A')
# Plot federal lands map
plot(vessel, add=TRUE,
     col = scales::alpha('#B4DEBC', 0.5), # RFF green
     border=NA)

# Add state lines again on top of federal lands, clear fill
plot(states_sp, add=TRUE,
     lwd=ifelse(use.png, 1.2, 1),
     col = 'transparent',
     border = '#687D8A')

# Plot points (equal size)
plot(points_spdf,
     pch = pch_lookup[points_spdf$fed.type],
     col = col_lookup[points_spdf$fed.type],
     cex = pch.cex_lookup[points_spdf$fed.type],
     add=TRUE)

# Add a box around the plot
box()

# Legend for colors
pt.reorder = c(4,3,2,1)
level.names = levels(points_spdf$fed.type)[pt.reorder]
level.names = gsub('Oil', 'Oil, ', level.names)
level.names = gsub('Gas', 'Gas, ', level.names)
legend(x=-127, y=ifelse(use.png, 26.5, 28),
       cex=1.25*ifelse(use.png, 1.5, 1),
       title='Wells', title.adj=0.2,
       legend = level.names,
       pt.cex = c(rep(ifelse(use.png, 3, 2), 2), rep(3*7.5/9, 2)),
       col = alpha(col_lookup[levels(points_spdf$fed.type)][pt.reorder], 0.6),
       pch = pch_lookup[levels(points_spdf$fed.type)][pt.reorder],
       bg='white')

legend(x=-127, y=ifelse(use.png, 29.5, 31),
       legend = 'Federal Lands',
       cex=1.25*ifelse(use.png, 1.5, 1),
       pt.cex = 2*ifelse(use.png, 1.5, 1),
       col = scales::alpha('#B4DEBC', 0.5), # RFF green
       pch = 15, bg='white')

dev.off()

difftime(Sys.time(), run.start)
rm(dt.map)
rm(points_spdf)
rm(points_sp)
rm(locs.map)
rm(vessel)
